\documentclass[12pt,UTF8]{ctexart}
\title{\textbf{最优化理论与算法课程第2次作业}}
\author{马哲辉}
\date{\today}
\usepackage[left=1.25in,right=1.25in,top=1in,bottom=1in]{geometry}

\usepackage{amsmath}
\usepackage{hyperref}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{ulem}
\usepackage{amssymb}
\usepackage{float}
\usepackage{caption}
\usepackage[dvipsnames]{xcolor} % 更全的色系
\usepackage{listings} % 排代码用的宏包



\begin{document}
\maketitle

\tableofcontents
\begin{center}
\includegraphics[scale=0.6]{题目.png}
\end{center}

\section{课程题目}
\subsection{最速下降法}

\begin{align*}
    f(x) & = x_{1}^{2}+2x_{2}^{2}-2x_{1}x_2\\
    \nabla f(x) &= \begin{bmatrix}
    2x_{1}-2x_{2}\\
    4x_{2}-2x_{1}
    \end{bmatrix}\\
    H &=\begin{bmatrix}
        \frac{\partial^2 f(x)}{\partial x_{1}^2} & \frac{\partial^2 f(x)}{\partial x_{1}\partial x_{2}}\\
        \frac{\partial^2 f(x)}{\partial x_{1}\partial x_{2}}&\frac{\partial^2 f(x)}{\partial x_{2}^2}
    \end{bmatrix}  = \begin{bmatrix}
        2 & -2\\
        -2 & 4 
    \end{bmatrix} \\
    \alpha_k & = \frac{(\nabla f(x))^T \nabla f(x)}{(\nabla f(x))^T H \nabla f(x)}
\end{align*} 

\lstset{
	language = matlab,
	backgroundcolor = \color{white!10}, % 背景色：淡黄
	basicstyle = \small\ttfamily, % 基本样式 + 小号字体
	rulesepcolor= \color{gray}, % 代码块边框颜色
	breaklines = true, % 代码过长则换行
	numbers = left, % 行号在左侧显示
	numberstyle = \small, % 行号字体
	keywordstyle = \color{blue}, % 关键字颜色
	commentstyle =\color{green!100}, % 注释颜色
	stringstyle = \color{red!100}, % 字符串颜色
	frame = lrtb, % 用(带影子效果)方框框住代码块
	showspaces = false, % 不显示空格
	columns = fixed, % 字间距固定
	%escapeinside={} % 特殊自定分隔符：
	morekeywords = {as}, % 自加新的关键字(必须前后都是空格)
	deletendkeywords = {compile} % 删除内定关键字；删除错误标记的关键字用deletekeywords删！
}
\begin{lstlisting}[caption=最速下降, language=matlab]
    clear
    clc
    close all
    %% 初始化
    n = 100;%迭代次数
    fx = @(x) x(1)^2 + 2*x(2)^2 -2*x(1)*x(2);%函数
    sx2 = @(x1) (2-x1.^2)./x1;%边界
    dfx = @(x) [2*x(1)-2*x(2),4*x(2)-2*x(1)];%梯度
    syms x1 x2;
    x = [x1,x2];
    f_sym = fx(x);
    dfx_1_1 = diff(diff(f_sym,x1),x1);
    dfx_1_2 = diff(diff(f_sym,x1),x2);
    dfx_2_2 = diff(diff(f_sym,x2),x2);
    
    H = double([dfx_1_1,dfx_1_2;dfx_1_2,dfx_2_2]);
    
    sample.x = [];
    sample.dx = [];
    sample.y = [];
    x_iteration = repmat(sample,n,1);
    %% 绘制计算域
    nx = 100;
    x1 = linspace(0,sqrt(2),nx);
    bx2 = sx2(x1);%边界
    x2 = x1;
    for i = 1:1:nx
        for j = 1:1:nx
            y(i,j) = fx([x1(i),x2(j)]);
        end
    end
    pcolor(x1,x2,y);
    shading interp;
    hold on 
    plot(x1,bx2,'r');
    
    %% 计算
    x_iteration(1).x = [1,1];%起始点
    x_iteration(1).y = fx(x_iteration(1).x);
    x_iteration(1).dx = dfx(x_iteration(1).x);
    
    iteration = [1];
    error = norm(x_iteration(1).dx);
    Y = x_iteration(1).y;
    %figure(2)
    %for i = 2:1:n
    i = 1;
    while (error(end) >= 1e-10) && (i <= n)
        i = i+1;
        % clf
        %梯度下降
        g = x_iteration(i-1).dx;
        g = g';
        alpha = (g'*g)/(g'*H*g);
        x_iteration(i).x = x_iteration(i-1).x - alpha*x_iteration(i-1).dx;
    
        %边界条件
        x_iteration(i).x(1) = min(x_iteration(i).x(1),sqrt(2)).*(x_iteration(i).x(1)>0);
        x_iteration(i).x(2) = (x_iteration(i).x(2) < sx2(x_iteration(i).x(1)))...
            .*(x_iteration(i).x(2) -sx2(x_iteration(i).x(1)))+sx2(x_iteration(i).x(1));
    
        x_iteration(i).y = fx(x_iteration(i).x);
    
        %求梯度
        x_iteration(i).dx = dfx(x_iteration(i).x);
    
        %绘图
        % iteration = [iteration,i];
        % error = [error,norm(x_iteration(i).dx)];
        % Y = [Y,x_iteration(i).y];
        % subplot(2,1,1)
        % semilogy(iteration,error);
        % subplot(2,1,2)
        % semilogy(iteration,Y);
        hold on
        x_all = [x_iteration.x];
        plot(x_all(1:2:end),x_all(2:2:end));
        pause(0.001);
    end
\end{lstlisting}

\begin{figure}[htbp] 
\begin{center}
    \includegraphics[scale = 0.5]{1.1.1.png}
    \caption{迭代过程}
\end{center}
\end{figure}

\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 0.75]{1.1.2.png}
        \caption{梯度值与函数值}
    \end{center}
    如上图所示，最优解为$x_1 = 0,x_2 = 0,f(x) = 0$。
\end{figure}

\subsection{共轭方向法}
\begin{align*}
    A &=\begin{bmatrix}
        2 & -2\\
        -2 & A 
    \end{bmatrix} \\
    s_0 &= -\nabla f(x_0)\\
    s_{k+1} &= -\nabla f(x_{k+1})+\beta_k s_k\\
    \beta_k &= \frac{\| \nabla f(x_{k+1}) \|^2}{\| \nabla f(x_{k})\|^2}\\
    \alpha_k &= -\frac{[\nabla f(x_k)]^{T}s_k}{2s^{T}_{k}As_k}
\end{align*} 

\begin{lstlisting}[caption=共轭方向, language=matlab]
    clear
    clc
    close all
    %% 初始化
    n = 1000;%迭代次数
    fx = @(x) x(1)^2 + 2*x(2)^2 -2*x(1)*x(2);%函数
    sx2 = @(x1) (2-x1^2)./x1;%边界
    dfx = @(x) [2*x(1)-2*x(2);4*x(2)-2*x(1)];%梯度
    
    H = [2,-2;-2,4];%Hessian矩阵
    
    sample.x = [];
    sample.dx = [];
    sample.y = [];
    x_iteration = repmat(sample,n,1);
    %% 计算
    x_iteration(1).x = [1;1];%起始点
    x_iteration(1).y = fx(x_iteration(1).x);
    x_iteration(1).dx = dfx(x_iteration(1).x);
    
    s = -x_iteration(1).dx;%初始的sk
    iteration = [1];
    error = norm(x_iteration(1).dx);
    Y = x_iteration(1).y;
    figure
    i = 1;
    while(error(end)>=1e-30 && i<=n)
    %for i = 2:1:n
        i = i+1;
        clf
        s0 = s;
        %梯度下降
        g = x_iteration(i-1).dx;
        alpha = -(g'*s)/(2*s'*H*s);
        x_iteration(i).x = x_iteration(i-1).x + alpha*s;
    
        %边界条件
        x_iteration(i).x(1) = min(x_iteration(i).x(1),sqrt(2)).*(x_iteration(i).x(1)>0);
        x_iteration(i).x(2) = ((x_iteration(i).x(2) < sx2(x_iteration(i).x(1)))...
            .*(x_iteration(i).x(2) -sx2(x_iteration(i).x(1)))+sx2(x_iteration(i).x(1))).*(x_iteration(i).x(2)>0);
    
        x_iteration(i).y = fx(x_iteration(i).x);
    
        %求梯度
        x_iteration(i).dx = dfx(x_iteration(i).x);
        beta = (norm(x_iteration(i).dx)^2)/(norm(x_iteration(i-1).dx)^2);
        s = -x_iteration(i).dx+beta*s0;
    
        %绘图
        iteration = [iteration,i];
        error = [error,norm(x_iteration(i).dx)];
        Y = [Y,x_iteration(i).y];
        subplot(2,1,1)
        semilogy(iteration,error);
        yticks(10.^(-30:3:0));
        subplot(2,1,2)
        semilogy(iteration,Y);
        yticks(10.^(-60:5:0));
        pause(0.001);
    end
\end{lstlisting}


\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{1.2.1.png}
        \caption{梯度值与函数值}
    \end{center}
    如上图所示，最优解为$x_1 = 0,x_2 = 0,f(x) = 0$。
\end{figure}


\subsection{罚函数方法}


\begin{align*}
    f(x) &= x_{1}^{2}+2x_{2}^{2}-2x_{1}x_{2}+ \rho \left( \frac{1}{2-x_{1}^{2}-x_1 x_2}
    -\frac{1}{x_1}-\frac{1}{x_2}\right)\\
\end{align*} 

初始罚因子为1，每过50次迭代，罚因子衰减为原先的$\frac{1}{10}$。

\begin{lstlisting}[caption=内罚函数法, language=matlab]
    clear
clc
close all
%% 初始化
n = 1000;%迭代次数
fx = @(x) x(1)^2 + 2*x(2)^2 -2*x(1)*x(2);%函数
sx2 = @(x1) (2-x1.^2)./x1;%边界
%dfx = @(x) [2*x(1)-2*x(2);4*x(2)-2*x(1)];%梯度
penalty = @(x,rho) rho*(1/(x(1)^2+x(1)*x(2)-2)-1/x(1)-1/x(2)); %罚函数
fx_pen = @(x,rho) fx(x)+penalty(x,rho);%原函数加上罚函数
syms x1 x2 rho;
x = [x1,x2];
f_sym = fx_pen(x,rho);
dfx_1_1 = diff(diff(f_sym,x1),x1);
dfx_1_2 = diff(diff(f_sym,x1),x2);
dfx_2_2 = diff(diff(f_sym,x2),x2);

dfx = [diff(f_sym,x1);diff(f_sym,x2)];%梯度
H = [dfx_1_1,dfx_1_2;dfx_1_2,dfx_2_2];
%此时的hessian矩阵并不是常数，需要具体运算
% example = subs(H,[[x1,x2],rho],[[1,1.1],1])
% ex_val = double(example);

sample.x = [];
sample.dx = [];
sample.y = [];
x_iteration = repmat(sample,n,1);
%% 绘制计算域
nx = 100;
x_1 = linspace(-sqrt(2),sqrt(2),nx);
bx2 = sx2(x_1);%边界
x_2 = x_1;
for i = 1:1:nx
    for j = 1:1:nx
        y(i,j) = fx([x_1(i),x_2(j)]);
    end
end
pcolor(x_1,x_2,y);
shading interp;
hold on 
plot(x_1,bx2,'r');
hold on
%% 计算
x_iteration(1).x = [1.1;1.1];%起始点在边界之外
% x_iteration(1).x = [0.9;0.9];
x_iteration(1).y = fx(x_iteration(1).x);
%x_iteration(1).dx = dfx(x_iteration(1).x);
iteration = [1];
Y = x_iteration(1).y;
%figure
i = 1;
r = 1;%初始罚因子
x_iteration(1).dx = double(subs(dfx,[[x1,x2],rho],[x_iteration(1).x',r]));
s = -x_iteration(1).dx;%初始的sk
error = norm(x_iteration(1).dx);
while(error(end)>=1e-30 && i<=n)
%for i = 2:1:n
    i = i+1;
    clf
    s0 = s;
    %求Hessian矩阵
    Hess = double(subs(H,[[x1,x2],rho],[x_iteration(i-1).x',r]));
    %梯度下降
    g = x_iteration(i-1).dx;
    alpha = -(g'*s)/(2*s'*Hess*s);
    %alpha = 0.001;
    x_iteration(i).x = x_iteration(i-1).x + alpha*s;

    x_iteration(i).y = fx(x_iteration(i).x);

    %求梯度
    x_iteration(i).dx = double(subs(dfx,[[x1,x2],rho],[[x_iteration(i).x'],r]));
    beta = (norm(x_iteration(i).dx)^2)/(norm(x_iteration(i-1).dx)^2);
    %s = -x_iteration(i).dx+beta*s0;
    s = -x_iteration(i).dx;

    %罚因子衰减
    if mod(i,50)==0
        r = r/10;
    end
end
\end{lstlisting}

\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{1.3.1.png}
        \caption{迭代过程}
    \end{center}
    如上图所示，从边界外出发，穿过了约束边界，达到最优解。最优解为$x_1 = 0,x_2 = 0,f(x) = 0$。
\end{figure}

\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{1.3.2.png}
        \caption{梯度与函数值}
    \end{center}
\end{figure}

\subsection{多目标方法}

\begin{align*}
    f_{1}(x) & = x_{1}^{2}+2x_{2}^{2}-2x_{1}x_2\\
    f_{2}(x) & = 2*(x_{1}-1)^2+(x_{2}-0.5)^2\\
    F(x) & = f_{1}(x)+f_{2}(x)
\end{align*}

\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{1.4.1.png}
        \caption{梯度与函数值}
    \end{center}
    最优解为$x_1 = 0.8125,x_2 = 0.4375,f(x) = 0.4062$。
\end{figure}

\section{生活题目}
小张家每天需要规划饮食和健身活动，以确保家庭成员的健康水平和控制家庭开支。他们有两个主要目标：一是提高家庭成员的健康水平，二是尽量减少不必要的开支。

\textbf{参数定义}

x：用于购买健康食品的预算（如新鲜蔬菜、水果、优质蛋白质等）

y：用于健身活动的预算（如健身房会员费、运动装备等）

\textbf{目标函数}

最大化健康水平$f_{1}(x_1,x_2) = -(x-100)^2-(y-150)^2+10000$

最小化不必要开支$f_{2}(x_1,x_2) = (x-50)^2+(y-100)^2$

\textbf{约束条件}

总预算限制$x+y \leq 300$

非负性$x \geq 0, y \geq 0$

\subsection{最速下降法}
\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{2.1.1.png}
        \caption{梯度与函数值}
    \end{center}
    最优解为$x_1 = 100,x_2 = 150,f(x) = -10000$。
\end{figure}

\subsection{罚函数方法}

\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{2.2.1.png}
        \caption{梯度与函数值}
    \end{center}
    最优解为$x_1 = 100,x_2 = 150,f(x) = -10000$。
\end{figure}

\subsection{多目标方法}

\begin{figure}[htbp] 
    \begin{center}
        \includegraphics[scale = 1]{2.3.1.png}
        \caption{梯度与函数值}
    \end{center}
    最优解为$x_1 = 75,x_2 = 125,f(x) = -7500$。
\end{figure}

\end{document}